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ABSTRACT 

The orbital evolution of a dust particle captured in a mean motion resonance with 
a planet in circular orbit under the action of the Poynting- Robertson effect, radial 
stellar wind and an interstellar gas flow of is investigated. The secular time derivative 
of Tisserand parameter is analytically derived for arbitrary orbit orientation. From the 
secular time derivative of Tisserand parameter a general relation between the secular 
time derivatives of eccentricity and inclination is obtained. In the planar case (the case 
when the initial dust particle position vector, initial dust particle velocity vector and 
interstellar gas velocity vector lie in the planet orbital plane) is possible to calculate 
directly the secular time derivative of eccentricity. 

Using numerical integration of equation of motion we confirmed our analytical 
results in the three-dimensional case and also in the planar case. Evolutions of ec- 
centricity of the dust particle captured in an exterior mean motion resonance under 
the action of the Poynting- Robertson effect, radial stellar wind for the cases with and 
without the interstellar gas flow are compared. 

Qualitative properties of the orbital evolution in the planar case are determined. 
Two main groups of the secular orbital evolutions exist. In the first group the eccen- 
tricity and argument of perihelion approach to some values. In the second group the 
eccentricity oscillates and argument of perihelion rapidly shifts. 
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1 INTRODUCTION 

For the analysis of dynamical evolution of dust particles in 
the vicinity of a star is necessary to take into account also 
non-gravitational effects. From the non-gravitational effects 
accelerations caused by the electromagnetic radiation and 
the corpuscular radiation (stellar wind) of the central star is 
most often considered. Influence of the electromagnetic radi- 
ation on the dynamics of dust grains is usually described us- 



ing Poynting- Robertson (PR) effect ( Poynting 1904 Robert- 
[son||ljff| |Klacka|p92| |2004| [Klacka et al.||2012b| )' For the 
corpuscular radiation is usually assumed radial propagation 
from the central star ( Whipple|1955 I. A more general form 
of the acceleration caused by the corpuscular radiation can 
be found in Klacka et al. (2012a I. When orbital periods of 



the particle and a planet are in a ratio of small natural 
numbers mean motion resonances (MMRs) can occur. If the 
dust particle is captured in an MMR, changes of the semi- 
major axis caused by small non-gravitational effects are bal- 
anced by the resonant interaction with the planet's gravity 



field. Jackson & Zook ( 1989 1 predicted a ring of dust par- 



ticles orbiting the Sun captured in mean motion resonances 
with planet Earth. This ring was confirmed by observations 
from satellit es IRAS (|Brownlee||1994| [Permott et al.||1994 | 
and C OBE ( [Reach et al.|1995| ). The paper [Jackson fc Zook| 
( 1989 1 was followed by many others who investigated be- 



haviour of dust particles captured in MMRs during the last 
three decades (e .g.lWeidenschilling fc Jackson'1993!'Marzari 



fc Vanzani|l994 iSidhchovsky & Ncsvorny 1994; Liou & Zook 



1995"; Liou, Zook fc 


Jackson 1995; Liou fc Zook 1997, 1999] 


Moro-Martfn fc Malhotra, 2002, Holmes ct al. 2003, Kuch- 


ner fc Holman|2003 


Deller fc Maddison|2005| Klacka, Koci- 


faj & Pastor||2005a 


b| Krivov et al.||2007| [Stark & Kuchner 


20081 Reach|2010||Mustill & Wyatt|2011||Ertel, Wolf & Rod- 



If the planet moves in a circular orbit around 
the Sun and mass of the particle is negligible in comparison 
with mass of the Sun and also in comparison with mass of 
the planet, we have a special gravitational problem of three 
bodies. This gravitational problem is called the circular re- 
stricted three-body problem (CR3BP) in celestial mechan- 
ics. An analytic expression for the secular time derivative of 
orbital eccentricity of the dust particle captured in an MMR 
in the planar CR3BP with the PR effect and the radial solar 



pavol.pastor@hvezdarenlevice.sk 



wind was found in Liou & Zook ( 19971. Due to the rotation 
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of the star stellar wind can be non-radial in general. Ac- 



cording to Helios 2 measurements ( Bruno et al. 2003 1 the 



angle between the radial direction and the direction of so- 
lar wind velocity is approximately constant (at least for the 
distances covered by observations) . The secular time deriva- 
tive of orbital eccentricity of the dust particle captured in 
an MMR in the planar CR3BP with the PR effect and such 



non-radial solar wind was calculated in Klacka et al. ( 2008 1 



and all possibilities of the secular eccentricity evolution were 
analytically determined in Pastor et al. ( 2009b I. The secular 



evolution of orbital eccentricity and argument of perihelion 
of the dust particle captured in an MMR in the planar cir- 
cular and elliptical restricted three-body problem with the 
PR effect was numerically investigated in [Pastor, Klacka &:| 
Komarl (|2009a|). 



Interstellar medium atoms penetrate into the Solar Sys- 
tem due to relative motion of the Solar System with respect 
to the interstellar medium. These approaching atoms form 
an interstellar gas flow (IGF). Influence of this IGF on mo- 
tion of dust particles orbiting a star was mentioned already 
in Whipple ( 1955 1. Motion of dust particles orbiting the Sun 



influenced by the IGF was analytically and numerically in- 



vestigated in Scherer (20001. When the IGF velocity vector 



lies in the orbital plane of the particle and the particle is 
under the action of the PR effect, radial solar wind, and 



an IGF, then the motion occurs in a plane. Scherer (20001 



has correctly described qualitative properties of the shift of 
perihelion in the planar case despite of the fact that his 
calculations contain several errors. Secular time derivatives 
of semimajor axis, eccentricity and argument of perihelion 



in the planar case were calculated in Klacka et al. ( 2009 1 



In Pastor, Klacka & Komar (20111 secular time derivatives 



of all Keplerian orbital elements for arbitrary orientation of 
the orbit with respect to interstellar gas velocity vector were 
calculated. The secular time derivatives of orbital elements 
were in [Pastor et al.| ( |2011[ ) derived under the assumptions 
(a) that the acceleration caused by the IGF is small com- 
pared to the gravitation of a central object, (b) that the 
speed of the IGF is large in comparison with the speed of 
the dust particle (speeds are determined with respect to the 
central object) and (c) that the speed of the IGF is large also 
in comparison with the mean thermal speed of the gas in the 
flow (Mach number 3> 1). Under these assumptions the IGF 
always causes decrease of the secular semimajor axis of the 
dust particle. This decrease of the secular semimajor axis 



was conflrmed analytically in Belyaev & Rafikov (20101 and 



using numerical integrations in Pastor et al. (2011 1, Marzari 
|fc Thebaul't] ( |2011| ) and |Marzari| ( |2012[ ). The acceleration of 
the dust particle caused by the IGF depends on the drag 
coefficient which is, for the given particle and the flow of in- 
terstellar gas, a specific function of the relative speed of the 



dust particle with respect to the interstellar gas (Baines 



Williams & Asebiomo 19651. The drag coefficient was in 
Scherer] (|2000|), [Klacka et al.| ( |2009l), [Pastor et al.| ( [20TT| , 
Marzari & Thebault[ (|2011[) and [Marzari[ ([2012[) taken into 



account as a constant. This is a consequence of assumption 
(c). The derivation of the secular time derivatives of all Kep- 



lerian orbital elements was generalised by adding variability 



of the drag coefficient during orbit in Pastor 



method applied in Belyaev & Rafikov (2010) for the secular 



(2012b I using a 



time derivative of semimajor axis. In view of assumptions 
this means that the secular time derivatives were derived 
under the assumptions (a), (b) and assumption (c) was re- 
placed by assumption that the mean thermal speed of the 



gas in the fiow is not close to zero (see Baines et al.[[1965 |. 
Also under these assumptions the secular semimajor axis 
of particle's orbit always decreases under the action of the 
IGF. The minimal and maximal values of the decrease of 
the semimajor axis were also determined in [Pasto r (2012b I. 
In this work we use the secular time derivatives of or- 



bital elements caused by an IGF derived in ^Pastor (|2012bl 



to obtain some basic properties of orbital evolution of the 
dust particle captured in an MMR with a planet in circular 
orbit. 



2 SECULAR ORBITAL EVOLUTION IN MMRS 

Within the framework of the CR3BP under the assumptions 
that the particle is far enough from the planet and mass of 
the planet is negligible with respect to the mass of the Sun F. 
F. Tisserand found a quantity which remains constant dur- 
ing the motion of the particle ( Tisserand|1896 
fc Clemcncc 1961| ) 



e.g. 



Brouwer 



(1) 



where a is the semimajor axis of the particle orbit, e the 
eccentricity of the particle orbit, i is the inclination of the 
particle orbit with respect to the planet orbital plane and 
ap is the semimajor axis of the planet orbit. If we take 
into account also non-gravitational effects, Tisserand param- 
eter will no longer be constant, in general. From the non- 
gravitational effects the electromagnetic radiation of the Sun 
in the form of the PR effect is most often considered. If we 
add Keplerian term of the PR effect to the central Keple- 
rian acceleration of the Sun, Tisserand parameter obtains 
the following form (cf. |Liou fc Zook|1997 l 



2afj 



(l-/3)a^(l 



cos 1/3 



(2) 



The parameter /3 is defined as the ratio of the electromag- 
netic radiation pressure force and the gravitational force be- 
tween the Sun and the particle at rest with respect to the 
Sun 



(3) 



Here, Lq is the solar luminosity, Q'p^ is the dimensionless 
efficiency factor for radiation pressure averaged over the so- 
lar spectrum and calculated for the radial direction (Qpr = 
1 for a perfectly absorbing sphere), c is the speed of light 
in vacuum, fi = GMq, G is the gravitational constant, Mq 
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is the mass of the Sun, and R is radius of the dust particle 
with the mass density g. Subscript /3 in Eq. ([2| denotes that 
osculating orbital elements are calculated using acceleration 
—fj-{l — I3)f/r'^ as the central Keplerian acceleration. Here, 
f is the position vector of the dust particle with respect to 
the Sun and r = |r|. 

The particle is in an MMR with the planet when the 
ratio of their mean motions is equal to the ratio of two small 
natural numbers. For an exterior q-th order MMR we have 
np/n = {p + q)/p and np/n = p/{p + q) for an interior g-th 
order MMR. Here, p and q are two small natural numbers, 
np is the mean motion of the planet and n is the mean 
motion of the particle. The special case mean motion 1/1 
resonance corresponds to q = 0. If the dust particle is cap- 
tured in the MMR, then the semimajor axis of the particle's 
orbit librates around a constant value. This can be written 
in the following form 



dap 
~dt 



= 



(4) 



Averaging in Eq. (|4| is over a period of the resonant li- 
bration. Using Kepler's Third Law we can obtain a relation 
between the semimajor axis of the planet and the semimajor 
axis of the particle captured in the MMR. We have 



G{Mq + Mp) 



2 3 
Upap 

2 3 

n ap , 



(5) 
(6) 



where Mp is the mass of the planet. Putting Eq. (|5| into 
Eq. ([6| (with assumption Mp <^ Mq) we find the relation 



afi = ap (1 - P) 



1/3 



np 
n 



2/3 



(7) 



Now, we admit that the dust particle is under the action 
of arbitrary non-gravitational effects for which secular time 
derivatives of the Keplerian orbital elements can be deter- 
mined. Because, the secular semimajor axis is constant, the 
gravitational influence of the planet (subscript G) must com- 
pensate changes of the semimajor axis caused by the non- 
gravitational effects (subscript EF) 



I dap 
\ dt 

We have also 



dap 



/ dap\ 
\dt L 



= 



dep \ 
dt I 

dip \ 



/dejlK 



1 

\ dt / j^p 



dip 
dt 



+ 



dip 
dt 



(8) 



(9) 



(10) 



For a particle with j3 = Q (e.g., an asteroid) in the planar 
CR3BP {ip = and dip /dt = 0) captured in an MMR the 
eccentricity oscillates around a constant value, because the 
Tisserand parameter (Eq. [T]) also oscillates around a con- 
stant value. Now we find a relation between the secular evo- 
lution of eccentricity and inclination for the dust particle 
under the action of the non-gravitational effects captured in 
an MMR. The total differential of the Tisserand parameter 



dt I 



dCr / dap \ ^ dCr / dep 
dap \ dt / dep \ dt 



+ 



dCr 
dip 



dip 
dt 



^ Ida^ 
a \ dt / pp 
I dea 



OCt f I dap 
dap W dt 

^dCrf/d^X ^ 

dep \\ dt / Q \ dt 

dCr ( I dip^ \ / dip^ 
dt / ^ \ dt 



+ 



(11) 



When the particle is far enough from the planet and the 
mass of the planet is negligible with respect to the mass of 
the Sun, then the planet gravitation does not change value 
of the Tisserand parameter. Hence 



dCr I dap 
dap \ dt I pp 

dCr /dAp\ 
dip \ dt / pp 



dCr /de^y 
dep \ dt / pp 



+ 



(12) 



Using the condition that the particle is in an MMR (Eq.|4| 
and Eq. ([12} in Eq. Jll]) we obtain 



+ 



dCr I dep 
dep \ dt 

dCr I dap \ 
dap \ dt I pp 

dCr I dip 
dip \ dt 



dCr 

dip 



dip \ 
dt I 



dCr 

dep 



dt 



+ 



(13) 



Eq. ( 13 1 represents the relation which must hold between the 



secular time derivatives of eccentricity and inclination of the 
dust particle captured in the MMR under the action of the 
non-gravitational effects. For the secular time derivative of 
eccentricity we obtain 



/de^\ _ dCr/dap I dap \ / dep \ 

\ dt / ' dCr/dep \ dt I pp ^\dt I pp 

_ dCr/dip I dip \ 
dCr/dep \dt / ^ ' 

Calculation of the partial derivatives yields 



(14) 



dep \ 
dt I 



■4 



2a pep 

' dap 
dt 



ap 



(1 - e|)i/2cosi^ \ap 



I dep 
\ dt 



3/2 



1 



4 



6/3 



tanip 



dip 
dt 



(15) 



Due to the derivation formulated for arbitrary non- 
gravitational effects with known secular time derivatives of 



orbital elements Eqs. ( 13 1 and ( 15 1 are generalisations of the 
results obtained by |Liou fc Zook| |T997l ) . 
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Figure 1. Evolutions of the orbit averaged semi-major axis, eccentricity, argument of perihelion, longitude of the ascending node, 
inclination and Tisserand parameter of a dust particle with R = 2 /im, g = I g/cm^, and Q'p^ = 1 captured in an exterior mean motion 
orbital 2/1 resonance with Neptune under the action of the PR effect, radial solar wind and IGF. 
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Figure 2. Evolutions of the orbit averaged time derivative of semi-major axis, eccentricity, argument of perihelion, longitude of the 
ascending node, incHnation and Tisserand parameter during the numerical solution of Eq. | |30| l depicted in Fig.[l](grey line). The evolution 
of the orbit averaged time derivative of Tisserand parameter is compared v^ith the secular time derivative of Tisserand parameter obtained 
from Eq. | |28| l (blacls line). The evolution of the orbit averaged time derivative of eccentricity is compared with the secular time derivative 
of eccentricity calculated from Eq. (pSl using the numerical values of {dip/dt) depicted in the bottom left panel (black line). 
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3 INFLUENCE OF AN IGF ON ORBITAL 
EVOLUTION OF DUST PARTICLES IN 
MMRS 

In this section we use the analytical theory from the previ- 
ous section in order to find equations for the secular time 
derivatives of the particle's orbital elements in the case when 
the particle, captured in an MMR, is under the action of the 
PR effect, radial solar wind and IGF. The secular changes 
caused by the IGF will be described using the secular time 
derivatives of orbital elements derived in |Pastor| ( [2012bJ . For 
the PR effect and radial solar wind the standard expressions 
will be used ( |Wyatt fc Whipple||1950[ [Klacka et al.||2012a l. 
The secular time derivatives of the orbital elements caused 
by these effects then are 



' da/} \ 
^ dt / 



' dep ' 
, dt , 



' du>f} \ 
K dt / pp 



'dflfj 
, dt 



2 + 3e^ 
%J a^(I-e2)3/2 



Pfi 



1 - 



1 + 



X + 



Sep 



c V ^Q^J 2a2(l-e2)i/2 



Vp 



Pfi 
m(1-/3) 



4 



E 



P/3 



^Vm(1-/3) 



3S> 

e/3 



+ C/3 



COS 

sin j;3 



X (Sp cos 0^/3 — Ja sinix)^) 
2smif3 Y/i(l-/3) 



3e^ sin uj/s aftgi 
1 - Vf 



(16) 



3//3 , f^,95.(/|-S|)(l-e|) 



(17) 



e^- 66^ + 4-4(1-6^; 
Se^sincj/j cr/jgi 



(18) 



dip \ 
dt I pp 



X (5^ cosaj;3 — sinix)^) 

CQi'^iVFCfj j Pfj 



(19) 



m(i-/?) 



3e/j cos Ljp apQi 

1-6? + 



"/3 



X (5/3 sintj^ + coscj^g) 



(20) 



where ujp is the argument of perihelion and is the longi- 



tude of the ascending node. The sums in Eqs. (16h(20l run 
over all particle species i in the IGF. r) is the ratio of solar 
wind energy to electromagnetic solar energy, both radiated 
per unit of time 



E' 



SW J 11 f'SW J *- 



(21) 



here, u is the speed of the solar wind with respect to the 
Sun, u — 450 km/s, nisw j and Usw j, are the masses and 
concentrations of the solar wind particles of j-th type at a 
distance r from the Sun, respectively. 77 = 0.38 for the Sun 
( [Klacka et al.|2012a[ ). coi is the drag coefficient for the dust 
particle at rest with respect to the Sun ( Baines et al.|1965 l 



1 



COi = 



f- 



+ 



+ ( 1+4 



+ (1 - S,) 



A)' 


-4, 








erf (so 










3soi 



(22) 



here, erf(si) is the error function erf(soi) = J^"' e~* 

dt, 5i is the fraction of impinging particles specularly re- 
flected at the surface (for the resting particles, there is 



assumed diffuse reflection) (Baines et al. 1965 Gustafson 



1994 1, Td is the temperature of the dust grain, and Ti is the 
temperature of the i-th gas component, soi is the molecular 
speed ratio for the dust particle at rest with respect to the 
Sun (Mach number) 



SOi = 



2kTi 



Vf , 



(23) 



rrii is the mass of the neutral atoms of type i, k is the Boltz- 
mann's constant, vf = {vfx, vfy,vfz) is the velocity of the 
IGF in the frame associated with the Sun, vf ~ \vf\, 7i is 
the collision parameter, 



7» 



: Ui—A 

m 



(24) 



here, rii is the concentration of the interstellar neutral atoms 
of type i, and A — -kFI? is the geometrical cross section of 



the dust particle with mass m, — apil — e 
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(25) 



^13 = , 

the parameter gi describes variability of the drag coefficient 
( Pastor||2012b[ ) 

/ 1 3 



9i ^ — 



1 

COi 



— Srii 



(26) 



and 

Si3 = (cosf2^ coscj/3 — sinf2/3 sina;;3 cosi^)«FX 

+ (sin coscj^g + cos fi^ sincj/j cosii3)vFY 

+ sincj^ sini/3iiFZ , 
IjS = (— cosil^j sincjg — sinf2;3 cosaj/3 cosi^)iiFx 

+ (— sin f2a sincj^ + cos fi/j cosoj^ cos i^)wFy 

+ coso;^ sini^Ui?^ 
C/3 = sinf2^ sini/jUFX 

+ cosipvpz ■ (27) 

If we use Eqs. ([l6|, ([l7| and ([20| in Eq. ([Ts}, we obtain 

/dCV\ ^ 9Ct /dep\ OCt ( dip\ 
\ dt / defi \ dt / difs \ dt / 



— cos Qis sin ipVpY 



2(l-e^)3/2 Iv) 




X (C/3 COSOJ^ sin j;3 — cosi/j) 

1 O" n Oi / , 

- 0-/3 cosi;3 + -^^( (7,3(5,3 smo;^ 

2 \ 

+ Ip coscj^) sin 1,3 — (5*^ + I^) cosia^ 



H 

a,3 



1 + 4 



(28) 



This is the relation between (de^/dt) and (di^/dt) which 
must hold during the orbital evolution of the dust particle 
captured in an MMR. When the position vector and the 
velocity vector of the dust particle lie in the planet orbital 
plane and the interstellar gas velocity vector lies also in this 
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Figure 3. Evolutions of eccentricity of a dust particle with R = 2 
fim, g = 1 g/cm^, and Q^^ = 1 captured in an exterior mean mo- 
tion orbital 2/1 resonance with Neptune under the action of the 
PR effect and radial solar wind for various initial eccentricities. 



plane, the dust particle's motion is coplanar. In this case we 
have 1/3=0 and (dip/dt) = 0. Using these assumptions we 
obtain for the secular time derivative of eccentricity from 
Eq. ([28| and Eq. Q 



/ 

\dt / 



1 + 



(1 



2 + 3e 



/3 



2(l-e2)3/2 np 



+ E 



COi-^iVp 



Pl3 



m(i-/?) 



3Il3_ 
2 VF 



+ 



1 + 



\9i_ 
2vl 



{Si +4 



(1 



4)V2 



6/9 



np 



1 + 



9i 



(29) 



4 NUMERICAL RESULTS 

In this section we want to use numerical solutions of equa- 
tion of motion in order to compare numerical results with 
the analytical results derived in the previous section and to 
find some main properties of dust particle orbital evolution 
in an MMR with a planet under the action of the solar ra- 
diation and IGF. 
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Figure 4. Evolutions of the orbit averaged semi-major axis, ec- 
centricity and argument of perihelion of a dust particle with R = 
2 /im, £1=1 g/cm^, and Q'p^ = 1 captured in an exterior mean 
motion orbital 2/1 resonance with Neptune under the action of 
the PR effect, radial solar wind and IGF in the planar case. 



4.1 Equation of motion 

If we consider the gravitation of the Sun, gravitation of one 
planet, PR effect, radial solar wind and IGF, then equa- 
tion of motion of the dust particle has form ( [Klacka et al.| 
2012a|b| [Baines et al.|1965| l 

dv 

dt ~ 



(1 - /3) efi 



5pr 



V ■ Er ^ V 

Br + - 

c c 



CDi^i\v — Vf\{v- Vf) 

i 

GMp 



\r - rp\ 
GMp 



Ar-rp) 



\rp\ 



-rp 



(30) 



where v = df/dt is the velocity of the dust particle, ep = f/r 



is the radial unit vector, cm is the drag coefficient (Baines 
et al.|1965| ) 



CDi 



1 1 



1 + 



erf(si) 



Si is the molecular speed ratio 



1 

4sf 



St = 



2kTi ' ' 



(31) 



(32) 



and fp is the position vector of the planet with respect to 
the Sun. 



4.2 



Comparing analytical and numerical results in 
the three-dimensional case 



From numerical solution of Eq. ( 30 1 the orbit averaged time 



derivative of Tisserand parameter can be obtained. We used 
the numerical solution of Eq. ( |30[ ) and Eq. ( |28[ ) for compari- 
son of the numerical and analytical values of the secular time 
derivatives of Tisserand parameter. In the IGF we consid- 
ered the primary and secondary populations of neutral hy- 
drogen atoms and neutral helium atoms. The primary popu- 
lation of neutral hydrogen atoms and neutral helium atoms 
represent the original atoms of the IGF which penetrate 
into the heliosphere. The secondary population of neutral 
hydrogen atoms comprises the former protons from the IGF 
that acquired electrons from interstellar H° between the bow 



2009 Alouani-Bibi 



shock and the heliopause (Frisch et al 
et al.|20TT |. We adopted the following parameters for these 
components in the IGF. ni = 0.059 cm"^ and Ti = 6100 



K for the primary population of neutral hydrogen ( Frisch et 



aL] |2009[ ), 



n2 = 0.059 cm"^ and T2 = 1 6500 K for the sec- 

ondary population of neutral hydrogen ( [Frisch et al.||2009l ) 
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Figure 5. Evolutions of the orbit averaged time derivatives of 
semi-major axis, eccentricity and argument of perihelion during 
the numerical solution of Eq. \30\ depicted in Fig. |4] (grey line). 
The evolution of the orbit averaged time derivative of eccentric- 
ity is compared with the secular time derivative of eccentricity 
obtained from Eq. (|29| (black line). 



and finally na = 0.015 cm~ and T3 = 6300 K for the neu- 
tral helium ( Lallement et al.]|2005 |. We have assumed that 
the interstellar gas velocity vector is equal for all compo- 
nents and identical to the velocity vector of the neutral he- 
lium entering the Solar system. The neutral helium enter 
the Solar system with a speed of about vf ~ 26.3 km s~^ 
( Lallement et al.|2005 l, and arrive from the direction of Xsd 
— 254.7° (heliocentric ecliptic longitude) and /3eci = 5.2° 
(heliocentric ecliptic latitude; Lallement et al. 20051. The 



components of the velocity in the ecliptic coordinates with 
the X-axis aligned towards the actual equinox vp ~ — 26.3 
km/s [cos(254.7°)cos(5.2°), sin(254.7°) cos(5.2°), sin(5.2°)] 
were transformed into a coordinate system with the xy-p\a.ne 
lying in the orbital plane of planet Neptune and the x-axis 
lying in the ecliptic plane. The orbital elements of the dust 
particle were determined in this reference coordinate system. 



The drag coefficients for Eq. ( 30 1 were calculated from Eq. 
(311. We assumed that the atoms are specularly reflected 
at the surface of the dust grain {Si = 1). The planet was 
initially located on the x-sods. Initial semimajor axis of the 
dust particle is, in general, computed from the relation i„ 
= a/3 nj,/n + A, where a/j ^ip/n is defined by Eq. ^ (with 
/3 calculated from Eq. [3| and A is a shift from the exact 
resonant semimajor axis. As the initial conditions for a dust 
particle with R = 2 fim, mass density g = I g/cm^, and 
Q'pr = 1, we used ap in = 2/1 + 0.001 AU, e,3 i„ = 0.5, 
1^/3 in = 60°, flji in = 40°, and in ~ 20°. The initial true 
anomaly of the dust particle was i„ = 345°. From the 
numerical solution of Eq. ( 30 1 we obtained evolutions of the 



orbit averaged orbital elements and the orbit averaged time 
derivatives of orbital elements. Values of Tisserand param- 
eter were determined in every time step from the oscular 
orbital elements. From these values evolution of the orbit 
averaged Tisserand parameter and the orbit averaged time 
derivative of Tisserand parameter was also obtained. Evolu- 
tions obtained from the numerical solution of Eq. ( |30[ l are 
depicted in Fig. [l] (black line) and Fig. [2] (grey line). We must 
note that the evolution duration 10^ years requires the size 
of an interstellar gas cloud 26.9 pc in the direction of the 
interstellar gas velocity vector (constant velocity with mag- 
rritude 26.3 km/s is assumed) and such situation cannot al- 
ways occur in the real galactic environment. To compare the 
numerical results with our analytical theory we calculated 
values of the secular time derivative of Tisserand parame- 
ter from Eq. ( |28| l using the numerically calculated secular 
orbital elements. The obtained evolution is depicted in the 
bottom right panel of Fig. |2] using black line. The lines are 
practically overlapping each other. We also calculated values 
of the secular time derivative of eccentricity from equation 



(see Eq. 28 I 



/de^\ 1 / dCr \ 

\ dt / dCr/deis \ dt / 



dCr/dip 
dCr/deis 



dip 
dt 



(33) 



where {dCr / dt) was calculated from Eq. ( |28[ ) and as {dip/ dt) 
we used the values obtained from numerical solution of Eq. 
( 30 1 and depicted in the bottom left panel of Fig. [2] Ob- 
tained evolution is shown in the top right panel of Fig. |2] 
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with black line. The analytical and numerical results are in 
excellent agreement. 



4.3 Comparing analytical and numerical results in 
the planar case 

In the planar case it is possible to compare directly the sec- 
ular time derivative of eccentricity obtained from numerical 
solution with the same result obtained from analytical the- 
ory. 

In the case when the dust particle captured in an MMR 
is only under the action of the PR effect and radial solar 
wind, secular time derivative of eccentricity is given by the 



first term in Eq. ( 29 1 , namely 



dep \ 
dt I 



1 + 



Qpr 



2(l-e2)3/2 np 



(34) 



IS an mcreasmg 



Function /(e) = (2 + 3e|) / [2(1 - e 
function of eccentricity. Function obtains values from /(O) — 
1 to lime^i l{e) = oo. Therefore, the secular eccentricity is 
always a decreasing function in an interior resonance. In an 
exterior resonance such value of eccentricity e/jc exists that 



l{ei3c) = 



"lie 



2(1 



np 
n 



p + q 



(35) 



For 6/3 < e^c the secular eccentricity asymptotically in- 
creases to e/sc and for e/s > epc the secular eccentricity 
asymptotically decreases to epc- In the special case mean 
motion 1/1 resonance secular eccentricity asymptotically de- 
creases to eisc = ( [Klacka et al.|2008| [Pastor et al.|2009b l. 
To visualise the secular evolution of eccentricity of the dust 
particle captured in an MMR only under the action of the 
PR effect and radial solar wind we numerically solved Eq. 
(301 without the third term. The planet was initially on 
the a;-axis. We used dust particle with R = 2 /im, £1=1 
g/cm'^ and Q'p^ = 1. The initial semimajor axis was ap i„ 
~ 0,0 2/1 + 0.001 AU. The initial eccentricities and argu- 
ments of perihelion were 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 
0.9 and 260°, 275°, 285°, 295°, 295°, 290°, 280°, 265°, 250°, 
respectively. The initial true anomaly was fp in — for all 
particles. Results are depicted in Fig. [3] The limiting value 
e/3c for mean motion 2/1 resonance obtained from Eq. (351 
is e/3c ~ 0.4812. As can be easily seen from Fig.[3]the secular 
eccentricity really approaches this value. 

An orbital evolution of the dust particle under the ac- 
tion of the PR effect, radial solar wind and IGF in the planar 
case is depicted in Figs. [4] and [5] We used the same param- 
eters for various gas components of the IGF, initial position 
of the planet and properties of the dust grain as in Figs, [l] 
and [2] The pi anar case was attained by a rotation of the 
interstellar gas velocity vector into the orbital plane of Nep- 
tune around axis lying in the orbital plane of Neptune and 
perpendicular to the interstellar gas velocity vector. The ro- 



tation angle was 3.7°. Initial orbital parameters of the dust 
grain were ap in = 2/1 + 0.001 AU, ep in = 0.3, ujp in 
— and ffj in = 160°. As can be seen in Fig. [4] in this 
case the secular evolution of eccentricity is not a monotonic 
function of time. The secular time derivatives of eccentric- 
ity obtained from this numerical integration and from Eq. 
( 29 \ are compared in Fig. [5] Again analytical theory and 
numerical solution are in excellent agreement. 



4.4 Qualitative properties of orbital evolution in 
the planar case 

In order to find qualitative properties of the dust particle 
orbital evolution in MMRs under the action of the PR ef- 
fect, radial solar wind and IGF numerical solution Eq. (301 



for many various initial conditions is needed. Calculation of 
the variable drag coefficients according to Eq. (31 1 is very 



time consuming therefore we used an approximation that 
the drag coefficients are constant. If the inequality \v\ <C vp 
holds during orbit, then this approximation is usable (see 
Pastor 2012b I. Many sets of initial conditions in Eq. (301 



with the constant drag coefficients were numerically solved 
in order to study evolution of the dust particle in an exte- 
rior MMR. One such set of numerical solutions is depicted 
in Fig. |6] We used the same properties of the IGF, initial 
position of the planet and properties of the dust grain as in 
Figs. |4] and [5] Initial orbital parameters of the dust grain 
were ap i„ = ap 2/1 + 0.001 AU, ep G {0.1, 0.2, 0.9}, 
(^p in £ {0, 5°, 10°, 355°} and fp in ~ 0. Therefore, we 
obtained 648 individual evolutions. Not all 648 evolutions 
are depicted in Fig. [6] An evolution is depicted in Fig. |6] 
only if the dust particle was captured and remained in the 
resonance longer than 10^ years. Each evolution depicted in 
Fig. [6] is represented by a marker in Fig. [7] Position of each 
marker corresponds to the initial eccentricity and initial ar- 
gument of perihelion. Colour of each marker corresponds to 
the capture time for given evolution. White colour would 
theoretically correspond to zero capture time and the dark- 
est shade of a given colour corresponds to the capture time 
2 X 10^ years. A place in Fig. [t] without marker corresponds 
to the initial eccentricity and initial argument of perihelion 
for which the dust particle was not captured in the reso- 
nance or remained in the resonance shorter than 10^ years. 
393 evolutions is depicted in Fig. [6] From many such sets of 
numerical solutions as set depicted in Fig. [6] we found that 
two main groups of orbital evolution in an exterior MMR 
under the action of the PR effect, radial solar wind and IGF 
exist. As was already mentioned all evolutions depicted in 
Fig. [6] have the initial true anomaly equal to zero. For sets 
with different initial value of the true anomaly are these two 
main groups preserved. These two groups have also longest 
capture times. 

Orbital evolutions in the first group are characterised by 
approach of the eccentricity and argument of perihelion to 
some "constant" values. We found that these values are not 
exactly constant but in comparison with the evolution be- 
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Figure 6. Evolutions of tlie orbit averaged Keplerian orbital elements of a dust particle with R = 2 /xm, g = 1 g/cm^, and Qp^ = 1 
captured in an exterior mean motion orbital 2/1 resonance with Neptune under the action of the PR effect, radial solar wind and IGF 
in the planar case. The initial values of orbital parameters were i„ = 2/1 + 0.001 AU, e^g S {0.1, 0.2, 0.9}, ui/j i„ G {0, 5°, 
10°, 355°} and i„ = 0. Depicted are only the evolutions with the capture time longer than 10^ years (see Fig.l7|l. 
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Figure 7. The initial eccentricities and arguments of perihelion of the evolutions depicted in Fig. |6] Colour of each marker, on a position 
given by the initial eccentricity and argument of perihelion, corresponds the capture time of obtained evolution. The capture times are 
from to 2 X 10® years and the colours are from white to the darkest shade of a given colour, respectively. If the capture time is shorter 
than 10^ years, then the corresponding position is empty. The initial conditions of evolutions belonging to the second group (evolutions 
with the fast monotonic shift of perihelion and oscillations of eccentricity) are marked with triangles (see text). 



fore this "stabilisation" are relatively slowly changing. The ment of perihelion ujfj ~ -57.16° in Fig. [6] Many evolutions 

stabilised direction of the position vector of the orbit's per- approach to this value of u/s as can be seen in the bottom 

ihelion is almost parallel with direction of the IGF velocity left panel of Figjm //3 = in the parallel case (see Eq. 27 1. If 

vector in Fig. [6] The parallel case corresponds to the argu- we use in Eqs. ( |16[ )-(|20| condition — 0, we obtain the sec- 
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ular time derivatives of orbital elements caused by the PR 
effect, radial solar wind and IGF described by a constant 
acceleration. If the dust particle is only under the action of 
the solar gravity (which can be reduced using the radial Ke- 
plerian term of the PR effect) and the constant acceleration 
caused by the IGF, then the secular orbital motion can be 
completely solved analytically (see Pastor|2012a Belyaev & 



|Rafikov 2010"). Perhaps, the following idea can appear: for 
constant acceleration caused by the IGF (tj^g = 0) the stabil- 
isation of eccentricity would correspond to trivial solution of 



Eq. (29l /;9 = and e/3 = e,3c (Eq. 351. Therefore, we solved 
Eq. ( 30 1 with acceleration caused by the IGF described by 



a constant vector (dependence on particle velocity in the 
acceleration was not included) for a set of initial conditions 
identical to the set used in Fig. [6] The stabilisation occurred 
also in this case. However, the argument of perihelion into 
which perihelia approached was not equal to -57.16° but 
-63°. Another important difference was in the evolution ec- 
centricity. The eccentricities approached to ~ 0.61 and 
not to e/3c ~ 0.4812. Therefore, at least in this simplified 
case, the stabilisation does not correspond to trivial solu- 



tion of Eq. (29l = and ep = epc- 

Orbital evolutions in the second group are characterised 
by a fast monotonic shift of perihelion and oscillations of 
eccentricity. Typical examples of this group are depicted 
in the bottom right panel of Fig. |6] as rapidly decreasing 
curves. The initial conditions which lead to this behaviour 
are marked with triangles in Fig. [7| 



5 CONCLUSION 

We investigated the orbital evolution of a dust particle cap- 
tured in an MMR with a planet in a circular orbit under the 
action of the PR effect, radial stellar wind and IGF. The sec- 
ular time derivative of Tisserand parameter is analytically 
derived for arbitrary orbit orientation using previously de- 
rived secular time derivatives of Keplerian orbital elements 
caused by the PR effect, radial stellar wind and IGF. From 
the secular time derivative of Tisserand parameter a rela- 
tion between the secular time derivatives of eccentricity and 
inclination can be obtained. In the planar case we derived 
directly the secular time derivative of eccentricity. 

We numerically solved equation of motion of the dust 
particle in order to compare the analytically derived re- 
sults with the numerically obtained results in the three- 
dimensional case and also in the planar case. In both cases 
analytical and numerical results are in excellent agreement. 
This implies that the theory developed in |Liou &i Zook| 
|1997^ for the PR effect and radial solar wind can be gen- 
eralised for arbitrary non-gravitational effects with known 
secular time derivatives. 

If the dust particle is captured in an exterior MMR 
in the planar CR3BP with the PR effect and radial stellar 
wind, then the secular eccentricity is a monotonic function 
of time. If we take into account also an IGF, complicated be- 
haviour with oscillations shows. This is most likely caused 



by dependence of the secular time derivative of eccentricity 
on the argument of pericentre due to directional character 
of the acceleration caused by the IGF. However, qualitative 
properties of the secular evolution of eccentricity and argu- 
ment of pericentre can be determined. We found that two 
main groups exist. In the first group the eccentricity and ar- 
gument of pericentre approach to some values. In the second 
group the eccentricity oscillates and argument of pericentre 
rapidly shifts. 
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